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Abstract 



We revisit the isocurvature density perturbations induced by quantum fluctu- 
ations of the axion field by extending a recently developed analytic method and 
approximations to a time-dependent scalar potential, which enables us to follow the 
evolution of the axion until it starts to oscillate. We find that, as the initial mis- 
alignment angle approaches the hilltop of the potential, the isocurvature perturba- 
tions become significantly enhanced, while the non-Gaussianity parameter increases 
slowly but surely. As a result, the isocurvature constraint on the inflation scale is 
tightened as H in i < 100 GeV for the axion decay constant f a < 10 10 GeV, near the 
smaller end of the axion dark matter window. We also derive useful formulae for 
the power spectrum and non-Gaussianity of the isocurvature perturbations. 
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1 Introduction 

The identity of dark matter is one of the central issues in cosmology and particle physics. 
The QCD axion is an interesting and plausible candidate for dark matter; it is a Nambu- 
Goldstone (NG) particle associated with the spontaneous break down of the Peccei-Quinn 
(PQ) symmetry introduced to solve the strong CP problem in QCD [TJ Ej. As the axion 
settles down at the potential minimum, the CP violating phase 9 is set to a vanishingly 
small value, solving the strong CP problem. 

The dynamical relaxation of the CP phase necessarily induces coherent oscillations of 
the axion, which contribute to dark matter, since the axion is stable in a cosmological 
time scale for the decay constant in the axion window [3], 

10 9 GeV <f a <10 12 GeV. (1) 

The lower bound comes from astrophysical constraints including the cooling argument of 
globular-cluster stars, and the upper bound from the requirement that the axion density 
should not exceed the observed dark matter density for the initial misalignment angle of 
order unity. Once the fine-tuning of the initial position is allowed, the upper bound can 
be relaxed to e.g. the GUT or string scale. 

One of the features of the axion is that its quantum fluctuation during inflation natu- 
rally induces an almost scale-invariant isocurvature density fluctuation, which would leave 
a distinctive imprint on the CMB spectrum. The observed CMB spectrum is known to 
be fitted extremely well by a nearly scale-invariant adiabatic density perturbation, and a 
mixture of the isocurvature perturbations is tightly constrained by the observations jU [5] . 
This constraint can be interpreted as an upper bound on the inflation energy density 
as the quantum fluctuation of the axion field is set by the Hubble parameter during 
inflation J6j [7]. 

The statistical information contained in the density fluctuations can be exploited by 
estimating higher order correlation functions. The non-Gaussianity of the isocurvature 
density fluctuations has been studied from both theoretical [H HI [101 [TTJ [12] and ob- 
servational [T3J HH US] point of view. In particular, the non-Gaussianity of the CDM 
isocurvature perturbation induced by the axion was estimated under an approximation 



assuming a quadratic potential for the axion [8]. This approximation however breaks 
down as the initial misalignment angle deviates from the CP symmetric vacuum, and 
an anharmonic effect must be taken into account. The anharmonic effect on the axion 
abundance and the isocurvature power spectrum was studied in Refs. [161 El HI EE1 EQ] . 
However, while there is a consensus on the anharmonic effect on the axion abundance, 
there are apparent disagreement as to how the isocurvature perturbations are affected by 
the anharmonic effect. 

Recently an analytic method and approximations to compute density perturbations 
induced by a curvaton field with a broad class of the potential was developed by Kawasaki 
and two of the present authors (Kobayashi and Takahashi) (KKT in the following) [2T] . 
Using this new method, we can analytically estimate the density perturbations generated 
in a curvaton mechanism 



[2H [25] for a potential significantly deviated from the 
simple quadratic one [211 EEJ E7] Jj The analytic and numerical results were found to agree 
with each other to a very good accuracy, typically, within several percent. This method 
was however limited to a time-independent potential for the curvaton. 

In this paper we extend the analytic method to a time- dependent potential, and then 
apply the method to the case of the QCD axion. We find it essential to define the onset 
of oscillations properly in order to evaluate the axion isocurvature perturbations. To this 
end we solve an attractor equation of motion for the axion until it starts to oscillate. 
This is the key to understand the anharmonic effect on the behavior of isocurvature 
perturbations. Based on the extended KKT method, we first correctly estimate the power 
spectrum and non-Gaussianity of the axion isocurvature perturbations for any axion initial 
position, resolving the aforementioned disagreement in the literature. We find that both 
power spectrum and non-Gaussianity are enhanced as the initial field value approaches 
the hilltop of the potential, thus giving an extremely tight constraint on the inflation 
scale for the axion decay constant f a = (9(10 9-10 ) GeV, near the smaller end of the axion 
dark matter window (pQ). We will also provide useful formulae for the power spectrum 



4 In the hilltop limit, the density fluctuations receive a dominant contribution from the fluctuations 
of the timing when the curvaton starts to oscillate, rather than from those of the oscillation amplitude, 
in contrast to what is usually assumed for the quadratic potential. To our knowledge, this observation 
was first made in Ref. [28], and its analytic understanding was developed in Ref. [27]. See also Ref. [29] 
for spectator field models in light of the spectral index after the Planck. 



and non-Gaussianity of the axion isocurvaure perturbations. 

Here let us summarize the comparison of our results with the works in the past. 
First of all, the axion abundance has been studied extensively in the past, including 
the hilltop limit, and we obtained results which are consistent with the previous works. 
Namely, the axion abundance is enhanced as the initial field value approaches the hilltop 
of the potential. However, concerning the power spectrum of the axion isocurvature 
perturbations, we find that there are some inconsistencies among the literature. As we 
shall see later, we find that the power spectrum of the axion isocurvature perturbation 
gets significantly enhanced toward the hilltop, and its behavior agrees well with the result 
of Ref. [17] valid near the hilltop. However, such an enhancement was not correctly taken 
into account in Ref. [20], which resulted in a much weaker constraint on the inflation 
scale for f a = 0(1O 9 ~ 10 ) GeV in their Fig. 1. In Ref. [30], it was pointed out that the 
power spectrum can be suppressed or enhanced by the anharmonic effects depending on 
the parameters. We do not confirm such suppression in our analysis. Furthermore, their 
enhancement factor is weaker than what we find. 

The rest of this paper is organized as follows. In Sec. [2] we summarize the basic 
properties of the QCD axion. In Sec. El we give a brief review of the recently developed 
analytic approach and extend it to include the time-dependent potential. We apply this 
analysis to the QCD axion, and calculate the derivative of the e-folding with respect to 
the axion field value at the horizon exit. In Sec. HI using the results in Sec. EJ we estimate 
the power spectrum and non-Gaussianity of the isocurvature perturbations, and derive 
constraints on the inflation scale as a function of the axion decay constant. We also provide 
useful formulae for the power spectrum and non-Gaussianity of the axion isocurvature 
perturbations; given the analytic expression for the axion dark matter abundance, one 
can easily estimate them using the formulae. The last section is devoted to conclusions. 

2 Axion Dark Matter 

We here briefly review the basics of the QCD axion dark matter. The axion is a NG boson 
associated with the spontaneous breaking of the PQ symmetry, which was introduced to 
solve the strong CP problem. Throughout this paper we assume that the PQ symmetry 



is already broken during inflation and is not restored after inflation_|, and that the PQ 



breaking scale remains unchanged during and after inflation. See Ref. [31] for the case 
where the PQ breaking scale significantly evolves after inflation. 

The axion has a flat potential protected by the PQ symmetry, which however is ex- 
plicitly broken by the QCD anomaly. Thus, as the comic temperature drops down to the 
QCD scale, Aqcd ~ 200 MeV, the QCD interactions become strong and non-perturbative, 
and the axion gradually acquires a finite mass from the QCD instanton effects. The axion 
potential is approximately given by 

V(a,T) = m a {Tffl (l - cos (£)) (2) 

where f a is the axion decay constant and m a {T) is the temperature dependent axion mass 
approximately given by [32] 

m a (T) « J Amo (t) f ° r T >y AqCD (3) 

[ m for T < A QCD 

with A « 0.1 and p ~ 4. We set the CP symmetric vacuum to be the origin of the axion. 
The axion mass at zero temperature is related to the decay constant as 

m « -2- — - ~ 6.0 x 10 b eV I , (4) 

1 + z J a \ Ja J 

where z = m u /rrid, rn n = 135 MeV, f n = 92 MeV, and we used the conventional value for 
z = 0.56 in the second equality. 

At a sufficiently high temperature, the axion mass is much smaller than the Hubble 
parameter. When the temperature becomes as low as T < 0(1) GeV, the axion starts to 
oscillate around the potential minimum. Throughout this paper we assume the radiation 
dominated Universe when the axion starts to oscillate. For a small initial misalignment 
angle, #* = a*// a <^ tt/2, the axion potential can be approximated with a quadratic 
potential. Then the axion abundance is given by 



5 The isocurvature perturbations are not generated if the PQ symmetry is restored during or after 
inflation. In this case, however, the topological defect such as the cosmic strings and domain walls are 
generated, and in particular, the PQ sector must be such that the domain wall number Now is unity. 



where h is the present-day Hubble parameter in units of 100 km s _1 Mpc -1 . One can 
see from this form that the observed dark matter abundance can be naturally explained 
by the coherent oscillations of the axion for the initial misalignment angle of order unity 
and for f a ~ 10 12 GeV. On the other hand, the initial misalignment angle must be finely 
tuned to be an extremely small value for f a ^> 10 12 GeV to avoid the overclosure of the 
Universe. 

If the axion field initially sits near the maximum of the potential, i.e., 6* ~ 7r, the 
commencement of the coherent oscillations can be significantly delayed. This is known 
as the anharmonic effect [161 E] • In the following we study its effect on the isocurvature 
fluctuations of the axion by using the recently developed analytic method. We shall see 
that the anharmonic effect leads to the enhancement of the abundance as well as the 
isocurvature perturbations of the axion. 

3 Analytical method of density perturbations 
3.1 Basic strategy 

Here we provide a basic strategy of the analytical method to compute the density pertur- 
bations generated by quantum fluctuations of a scalar field. See Ref. [21] for details. 

The density perturbations induced by a light scalar field depend on the scalar evolu- 
tion during and after inflation. Suppose that the scalar potential is approximated by a 
quadratic potential, and that its mass is sufficiently light during and for some time after 
inflation. Then the scalar field hardly evolves and stays more or less at the initial position 
until it starts to oscillate. When the Hubble parameter becomes comparable to the mass, 
it starts to oscillate, and importantly, the timing does not depend on the position. This is 
no longer the case for a potential of the general form, and one needs to follow the evolution 
of the scalar field until the commencement of coherent oscillations in order to compute the 
density perturbations. To this end, we first note that the scalar evolution can be actually 
well described by an attractor equation of motion. Then, by using the attractor equation 
of motion, we express the dependence of the e-folding number on the initial position of 
the scalar field. Finally we compute the power spectrum and non-Gaussianity parameter 
of the isocurvature perturbations, making use of the 5N- formalism [331 EH EHJ [36] . The 



most important ingredients are dN/da* and d 2 N/da 2 , where N is the e- folding number 
between the horizon exit and some time after the matter-radiation equality, and a* is the 
field value at the horizon exit of the CMB scales. The purpose of the rest of this section 
is to express dN/da* and d 2 N/da 2 in terms of the scalar potential and the axion field. 

In Ref. [21], the scalar potential was assumed to be time-independent. Here we shall 
extend the attractor equation of motion to include a possible time-dependence of the 
scalar potential, which is an essential feature of the axion. 

3.2 The attractor equation of motion 

Here we derive an attractor equation of motion of the following form, 

cHa + V'{a,t) = 0, (6) 

where c is a positive constant, a dot represents a ^-derivative, and a prime a partial 
differentiation with respect to the scalar field a. In the following discussion we identify 
a with the QCD axion, but most of the results in this section can be straightforwardly 
applied to a generic scalar field with a potential V(a, t) with explicit time dependence. 

The attractor equation for a time-independent potential was derived in Ref. [37] and 
also in Appendix A of Ref. [21]. The point is that, the scalar evolution can be well 
described by a first order differential equation when the curvature of the potential is 
much smaller than the Hubble parameter. The coefficient c is determined so that the 
attractor equation is consistent with the true equation of motion, a + 3H a + V = 0. 

To simplify our analysis we consider a potential of the form, V(a,t) = f(t)v(a), with 
f(t) = —xHf(t), where a; is a constant. Then the constant c in (J2J) should satisfy 

This shows that in a Universe with constant H/H 2 (i.e. constant equation of state 
parameter w = p/p), then as long as that the potential curvature is as small as \fv"\ ^C 
c 2 H 2 , the constant c is given by 

H 

c^3- — -x. (8) 



One can further check that the approximation (jBJ) with (JSj) is actually a stable attractor 
for c > 0. 

In the case of axion, the temperature dependence is given bjc 

f(t) oc T~*, (9) 

for T ^> Aqcd- This leads to 

f - = -2p^ = 2pH, (10) 

i.e., x = —2p rj —8, where we have assumed radiation-domination, H = l/2t, in the 
second equality. (Having a constant x = —f/Hf greatly simplifies the analysis for axions, 
as we will soon see.) The coefficient c can be approximated with 

c ^ 5 + 2p, (11) 

where we used H = —2H 2 . The attractor equation of motion will be valid as long as 

< c 2 ~ 170. (12) 



fv" 



H 2 

Roughly speaking, the attractor equation of motion holds until the curvature of the po- 
tential becomes about 10 times as large as the Hubble parameter. 

3.3 The onset of oscillations 

The early stage of the scalar evolution can be described well by the attractor equation, 
which however breaks down at a certain point, and the scalar field starts to oscillate 
around the potential minimum. The purpose of this subsection is to define the timing of 
the commencement of oscillations, t = t osc , or equivalently, H = H osc , and to calculate its 
dependence on the initial position a*. Specifically, we will calculate d\nH 2 sc /da osc and 
da osc /da*. 

Setting the potential minimum around which the scalar oscillates to a = 0, the oscil- 
lations are considered to start when 

= «, (13) 



Ha 



6 We have confirmed that such approximation is indeed valid until the commencement of oscillations 
for the parameter region of our interest. 



where k is a constant of order unity. We will set it to be unity in the numerical calculations, 
but we have confirmed that our results remain almost intact as long as k is of order unity. 
Combined with the attractor equation, we obtain 

1 V 
HL = k^ ' (14) 

where it is assumed that the potential is an increasing function of a from the origin to 
the field values of interest. Differentiating both sides with respect to a osc , we obtain 

<} ■ 1- HL = -J- + -L (v'L + xfjp.^- 1 ( :i.r, ) 



da osc a v' \ osc sc f da 

•^""OSC "'OSC ^osc \ Jose UUj 

where / osc = /(^osc) and f osc = ^(aosc), an d ft should be noted that the potential explicitly 
depends on the time, which also depends on a osc . (In other words, a osc and t osc are related 
through Eq. ffT3|) .) Using 

^ = - — #ln^ 2 , (16) 

da 4H da K ' 

which holds in the radiation dominated era, we obtain 



"• 



da, 
with 



lniJ 2 sc = — X(a osc ), (17) 



X{a osc ) = (l 


' 1) \v[ 


later use we rewrite Eq. (TT7|) as 




df 

>- /( 'OSC 


X{a osc ) 


OQ^osc 


■" osc "osc 



// 

osc 

/ 
osc 



(18) 



(19) 

If the scalar potential is quadratic, the function X vanishes. On the other hand, in the 
hilltop limit, X becomes much larger than unity. Thus, X is considered to represent the 
effect of the deviation of the scalar potential from quadratic one. 

Next we calculate da^c/da*. Integrating the attractor equation over a = a* ~ a 
and H = H * ~ H mr , we obtain 



osc 



l OSC5 






1 f HoBC f(H) 

-da = const. - / J -±-^-dH ) (20) 



v'(a) J cHH 



where f(H) actually means f(t(H)), and terms that are independent of a* are denoted 
as const. We differentiate both sides with respect to a* to obtain 

— — = (1 - KX(a osc )) . (21) 

da* v'{a*) 

The main results of this section are ( 1T71) and ( 12~TT) . which will be used later to calculate 

ON /da*. 

3.4 Axion number density 

Next we estimate the axion number density. Since the axion mass increases after the 
onset of oscillations, it is the number density that determines the relic abundance of the 
axion; the number density decreases in proportional to the inverse of the volume, and its 
ratio to the entropy density is fixed afterwards. The axion energy density at a later time 
can be estimated by multiplying the number density with the zero-temperature mass, m . 
Denoting the mass at the origin on the onset of oscillations as m osc , the number density 
is related to the potential and mass as 

^a,osc — • K^"- 1 ) 

^osc 

Precisely speaking, the number density is also dependent on the kinetic energy. However, 
one can see that the kinetic energy contribution becomes smaller as one goes to the hilltop 
region, since the onset of oscillations are delayed. This can be seen as follows. The kinetic 
energy can be estimated as 

K.E. = l -a 2 osc ~ y# 2 sc aL, (23) 

where (Tl3|) is used. On the other hand, the potential energy is 

P.E. ~ m 2 osc a 2 osc . (24) 

Thus, for the delayed onset of oscillations, H osc < m osc , the kinetic energy is smaller than 
the potential energy. 

Differentiating the number density with respect to a osc , one obtains 

lnn„ 



da 



a,osc 
osc 







(25) 



10 



where it should be noted that m osc has no explicit dependence on a osc , and so, m' osc = 0. 
Using m osc oc / osc and Eq. f lT9|) . one arrives at 

lnn ajOSC ~ _°S£ + __i L. (26) 



— ■ -u,uyc — i q 
OCl osc ^osc ^ Oi 

3.5 The e-folding number 

Now we are ready to calculate the derivative of the e-folding number with respect to a*, 
which is directly related to the primordial curvature perturbations in the SN formalism. 
We are interested in the e-folding number between the horizon exit of the cosmological 
scales and some time after the matter-radiation equality, and these times are represented 
by t* and t ea d- Specifically we take the slicings at t — t osc and t = t en d as the flat slicing 
and uniform density slicing, respectively. The e-folding number between t* and t en d is 
given by 

N = J H(t')dt' } (27) 

Let us split the integral into two pieces; 

ftosc 

N a = H(t')dt', (28) 

Ju 

/*^end 

Np = H{t')dt'. (29) 

It is useful to use the radiation energy density, or the entropy density, instead of the time. 
We adopt the entropy density, s, which scales as the inverse of the volume, even if the 
relativistic degrees of freedom changes with time; 

s + 3Hs = 0. (30) 

The e-folding numbers can be rewritten as 



(31) 



N a = 


1, 
= -- In s osc + const 


N p = 


J- , "Send 

= — In , 

*J "Snsc 
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where s osc = s{t osc ) and s en d = s(t en d), and the constant term in N a does not depend on 
a*. The radiation energy density p r is related to the entropy density as 



7Tsu 4 7c 2 g* ( 45 V -i 



4 



4 



We are interested in the derivative of the e- folding number with respect to a*. The varia- 
tion of a*, 5a*, is set by the Hubble parameter during inflation, H*, and the corresponding 
fractional change of the temperature, 5T osc /T osc , is much smaller than unity for the pa- 
rameters of our interest. Thus, we can practically neglect the change of g*(T osc ) and 
g* s (T osc ) due to the variation of a*. Then we obtain 

da* ~ 4da/ npr ' osc - 4da, osc ' [SS) 

9a, " 4 9a, ^ p nosc " 4 9a, ln Pr ' end + 4 da* ^ ^° sc ' ^ 

where we have assumed that the Universe was radiation-dominated and the axion and 
the other CDM components were negligible at t = t osc . The advantage of using the 
entropy density is that it becomes clear that the above formulae are still valid even if the 
relativistic degrees of freedom changes between t osc and tend- 
Now we estimate the first term in Eq. (134"]) . which turns out to be rather involved. To 
this end we define 



R 



Pc 

Pr 
Pa 
Pc 



(35) 

end 



(36) 



end 



where the total DM density is given by p c = p m + p a , and p m denotes the CDM component 
other than the QCD axion. Using R and r, we have (note that H en d on the final uniform 
density slicing is a constant that is chosen independently of a*) 
d d 

l n Pr,end = q hi (3H cnd M p — p m , ea d ~ Pa,end) 



da* da* 

I "Pm, end CPa,end 



Pr,end \ da* da* 

dR{\ -r) d , _ 3N v 

Rr —— ln im n a>osc e ") 

da* da* 

(37) 



-R(l - r)— lnp rie nd tj Rr d~ ln ( m ° na > osce 3Np ) ■ 



12 



We suppose that the entropy perturbation between p m and p r is not generated by the 
fluctuation of a*. (Hence we consider the CDM components other than the axion are 
produced while the axion has little effect on the expansion of the Universe.) Thus, 



9 ■lnf^l = o. cm 



namely, 



Od-t V S / end 



d Id 

—— lni?(l-r) = -7-5— Intend- (39) 

da* 4 da* 



So, we obtain 



d , ARr fd\nn a , osc n dNp 



ln Pnend = --r-^rni t ; ' -3^ (40) 



da* ' ' 4 + 3i?(l — r) \ da* da, 

Substituting this into Eq. (134"|) . we arrive at 



^iV/3 Rr dlnry^ 4 + 3g(l - r) d 2 

da* A + 3R da, 4(4 + 3/2) da* n osc ' l J 

To summarize, 

ON Rr d / 3 2 \ 

da* 4 + 3i?da* V 4 osc J ' 



X )-lIosc Iosc_ 3 __ (42) 

4 + 3iT ; < V^osc V 2 J a os J' V ; 

where we used Eq. (117)1 and Eq. (|26|) in the last equality, and X is understood as X(a osc ). 
This result is reasonable, since it vanishes in the limit of r — > 0, i.e., when the axion has 
negligible energy density. 

In order to estimate the non-Gaussianity of the isocurvature perturbations, we need 
to calculate d 2 N/da\. The following formulae are useful for this purpose; 

da, da, ' 

dRr dN 

_ = (4 + fl (3 + r) )^, (44, 

which can be shown by direct calculation. After long calculation, we obtain the second 
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derivative of N as 

d 2 N 16 + 8R(3 -r) + 9(1 - r)R 2 fdN\ 2 Rr d (A + 3R8N 



Oa'i {A + 3R)Rr \daj + 4 + 3Rda*\ Rr da J ,! " 

16 + 8R(3 - r) + 9(1 - r)R 2 fdN^ 2 



(A + 3R)Rr \da 



Rr (1 - kX)- 1 % 



A + 3R 



kX ' V<*c f V 'osc _ ( 3 _ X\ X 



1 - nXy u£ V*W V 2 



+ (1 - -X)- 1 ('% - ^ (1 - kX)\ (<± - (3 " X 



K v'/ J \v osc V 2/a. 

v-l "^osc /"^osc "^osc f~ X\ f X X 



(1 - kX)" 1 I2E I25S - I^2_ - 3 - - 



^* \ ^osc ^osc V "osc ^osc 

(46) 
Now we are ready to estimate the isocurvature perturbations and its non-Gaussianity. 

4 The axion isocurvature perturbations and its non- 
Gaussianity 

4.1 Power spectrum and bi-spectrum 

Let us start by discussing how the two-point and three-point correlation functions of the 
isocurvature perturbations can be expressed in the SN formalism (331 |3U ESI [36] , following 
Ref. [8]. 

The isocurvature perturbation S(x) is defined as 

S(x) = 3(Cc(aO-Cr(£)), (47) 

where C c (r) is the curvature perturbation on the slicing where the dark matter (radiation) 
is spatially homogeneous. Using the 5N formalism, £ c(r ) is obtained as the fluctuations 
in the number of e-folds between an initial flat slicing (which we take as the time of 
CMB scale horizon exit) and a final uniform-p c(r ) slicing, among different patches of the 
Universe. Since we are interested in the difference between ( c and ( r , we only need to 
compute the contribution to ( c arising from the axion field fluctuations 5a*. Thus (|4"T|) 
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can be expressed in terms of the derivative of the e-folding number as [8 

S(f) = 3 (^ a *( f ) + ^^(^( f )-(^( f ))) + •••)• ( 48 ) 

By taking the final uniform-p c slicing for iV to be deep in the matter dominated era, i.e., 
R — > oo, the slicing can be identified with a uniform total density slicing and thus we can 
use the results from the previous section. 

Here, one may wonder about the validity of using the SN formalism since the CDM and 
radiation components exchange energy with each other, as is clearly seen as the axion's 
temperature dependent mass. We remark that the energy density of CDM is negligibly 
tiny compared to the total density while the axion mass depends on temperature (i.e. 
T ^> Aqcd)- We suppose this to be the case also for the other CDM components having 
interactions with radiation, say, during their thermal production. Therefore the axion 
field fluctuations have negligible effect on the expansion of the Universe while the CDM 
and radiation interact with each other, allowing us to freely choose the initial flat slicing 
to be not only after the interactions are turned off, but also before then. 

The power spectrum and the bispectrum are usually defined in the momentum space 
as 

(S(h)S(k 2 )) = (27i) 3 P s (k 1 )5 3 (k 1 + k 2 ), (49) 

(S^S^Siks)) = (27r) 3 B s {k 1 , k 2 , k 3 )S 3 (h + k 2 + k 3 ), (50) 

where the isocurvature perturbation in the momentum space is defined by 

S{k) = [d 3 xe- il *S(x). (51) 

The momentum dependence of the bispectrum can be approximately factored out; 

B s (k h k 2 , fc 3 ) = /,J L SO) [Ps(ki)P s {k 2 ) + 2 perms] , (52) 

where f^ L is the non-linearity parameter for the isocurvature perturbationsjj 

Here we do not consider fluctuations of CDM components other than the axion, and further suppose 
that there are no mixing terms such as d 2 N/da*d(f>* where 4> is the inflaton field. 
8 /nl i s identical to fs used in Ref. |S]. 
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Similarly we define the power spectrum of the axion field fluctuations; 

(<M*i)M£)) = (2rr) 3 P 5a Ah)S(ki + h) (53) 

where 

PsaM = g- (54) 

The power spectrum and the non-linearity parameter of the isocurvature perturbation 
can be expressed as, to leading order in the field fluctuationsjj 



Ps(k) = 9(^-j P Sa M: 



where we consider the axion to have Gaussian field fluctuations. Let us repeat that in 
the above expressions, we substitute the results from the previous section deep in the 
matter-dominated era, i.e., R — >■ oo. 

The constraint from the Planck and the WMAP large-scale polarization data reads [5] 

a < 0.041 (95%C.L.), (56) 

where a represents the relative amplitude of the power spectrum of the isocurvature 
perturbation to that of the adiabatic one^ 

" = now (57) 

with k Q = 0.05 Mpc -1 . The power spectrum of the curvature perturbations are given 

by m 

A 2 ( (k ) = ^Pc(ko) * 2.2 xlO" 9 . (58) 



9 It is straightforward to take into account of higher-order contributions, which can be relevant if 5a* 
is not sufficiently small compared to a*. This is the case if one considers a large value of -ff; n f and/or 
a small value of r. Note however that the expansion of SN in terms of 5a breaks down for 5a„ ~ a*. 
When the potential can be approximated by a quadratic one, an order-of-magnitude estimate of P$ is 
still possible by simply replacing 9% with 91 + <jg, where ag = H* / (2ir f a ) ■ 
10 Note that a defined here corresponds to /3i so /(l — /3iso) in Ref. [5]. 
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The bound on the non-linearity parameter f^ L was derived using the WMAP 7-yr 
data 



a 2 fS?° ] = 40 ± 66, 



(59) 



for the scale-invariant isocurvature perturbations, assuming the absence of non-Gaussianity 
in the adiabatic perturbations 







4.2 Axionic isocurvature perturbations and its non-Gaussianity 

Now we are ready to obtain the analytical expression for the isocurvature perturbation 
and its non-Gaussianity generated by the axion, by combining the results in the preceding 
sections. 

Let us first calculate the power spectrum of the isocurvature perturbation. To simplify 
the expression, we define A|(fc) = ^-P^A;). Using Eq. ( H2|) . we then obtain 

A s (/c) = 



1 - kX vL 



V -^-(l- X - 



X 



2^ 



(60) 



where we have taken the limit of R — y oo. In the hilltop limit, X becomes much larger 
than unity, and both v' osc and v[ approach 0. Actually, however, v' osc decreases much 
more slowly than v'^ [21J. Thus, the power spectrum of the isocurvature perturbation gets 
enhanced significantly by the factor v' osc /vl , in the hilltop limit. 
The non-Gaussianity parameter /^ L is given by 



f\ 



(ISO) 
NL 



— r 1 
+ - 



-3 



+ 



^OSC [q_ ^ 



X 



x 

2/ a, 
-l 



X 



-1 r 



V. 



kX' 

1-kX 

2 



+ 



(',; 



v' 



s "f 



v 
X' 



;i - kX) 



x 



.(61) 



In contrast to the power spectrum, there is no huge enhancement due to 1/v'^ in the hilltop 
limit; the increase of /^ L is much milder and is mainly due to l/fo SC - Note also that 
/nl does not depend on the inflation scale, while the power spectrum As does. 

In order to see the consistency of the above result with the previously known results, 
let us consider a quadratic potential, for which X = and X' = 0. This corresponds to 



We expect that the constraint on f^ L can be improved considerably by using the Planck data. 
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the case where the initial misalignment angle is small. Then one can easily check 

/r» ->£-!, (62) 

in the limit of the quadratic potential, v(a) oc a 2 . This coincides with the result, f^ L ~ 
l/2r for r < 1 in Ref. [8]. 

Using the above relations (160]) and (161]) . let us estimate the isocurvature perturbations 
and its non-Gaussianity of the axion dark matter. To this end, we need to solve Eq. (1T5]) 
(or Eq. (120]) ) in order to evaluate the commencement of oscillations, which determines a osc , 
v osc as well as X(a osc ). We have numerically solved Eq. f fT3|) for a osc , and then substitute 
it to the above relations to obtain As and f^ L . We present these semi-analytic results 
in the following. 

In Fig. [H we show the contours of the fraction of the axion to the total dark matter 
density, r = p a /(p a + Pm) (see Eq. ({55]) ). and the non-linearity parameter /^ L , in 
the plane of the initial misalignment angle a*// a and the axion decay constant f a - In 
the region above the top (red) solid line, the axion abundance exceeds the observed dark 
matter abundance. We can see that, for a given value of f a , the axion abundance increases 
as the initial position approaches the hilltop, a*/ f a = it. Accordingly, the observed dark 
matter abundance is realized at a smaller value of f a < 10 10 GeV for the hilltop initial 
condition, 1 — a*/n f a < 10~ 5 . The contours of f^ L show that they are parallel to those 
of r for a*// < 1. This is because, as we have seen above, f^ L is determined simply by r 
when the scalar potential can be approximated with the quadratic potential. On the other 
hand, as the initial position approaches the hilltop, the value of f^ L mildly increases 
along the contours of r. This mild increase is due to the delayed onset of oscillations, 
which is one of the features of the hilltop curvaton [21] . 

The inflationary scale, i^mf, is bounded from above by the isocurvature constraint 
(1551) . The contours of the upper bound on log 10 (if in f/ GeV) are shown in Fig. [2J Since 
the isocurvature perturbation is significantly enhanced towards the hilltop, the constraint 
on if in f becomes tight. For instance, for 1 — o*/7r/ a < 10~ 8 , H lTL t must be smaller than 
~ 1 GeVo On the other hand, the non-Gaussianity constraint ( 15 9 p gives much weaker 



12 We have checked that the axion does not fall into wrong vacuum by climbing over the hilltop of the 
potential due to the quantum fluctuations, as long as the isocurvature constraint on H in f is satisfied. 
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Figure 1: The contours of the fraction of the axion density to the total dark matter 
density (thick solid lines), r, and the non-linearity parameter (dashed lines) f^ L . The 
contours of r correspond to r = 1,0.1, 0.01, 10 -3 and 10~ 4 , from top to bottom. (The 
contours of r = 0.01, 10 -3 , 10~ 4 are not shown in the right panel.) 



constraints in the hilltop region. One can see this by noting that f^ L is at most ~ 100 
in the hilltop when r > 0.1, leading to c^/nl < (9(0.1). If the axion is a tiny fraction 
of dark matter (say, r < 10~ 4 ), the non-Gaussianity constraint becomes important, as 
shown in Ref. [8]. 

Lastly, we show the constraint on H in f as a function of the decay constant f a in Fig. |3j 
For comparison, we show two constraints, a < 0.041, and 0.001. The first one represents 
the current constraints from the Planck and WMAP large-scale polarization data, while 
the latter one is just for showing how the constraint changes as the constraint on a 
improves. We can see that the constraint becomes extremely tight for f a < 10 10 GeV due 
to the anharmonic effect near the hilltop. The effect becomes milder as r decreases since 
the initial misalignment angles deviates from the hilltop. 



4.3 Useful expressions for As and /, 



(ISO) 

NL 



In order to see if the above semi-analytic estimate correctly describes the power spectrum 
and non-Gaussianity of the axion isocurvature perturbations, let us compare them with 
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Figure 2: The contours of the upper bound on the inflation scale, Log 10 (if; n f/GeV), 
from the isocurvature perturbation constraint ( |56|) (thick solid lines). The contours of r 
in Fig. []] are shown together. 



another (more direct) calculation. To this end, it is important to note that (i) the ac- 
tual abundance of the QCD axion is much smaller than that of radiation, which makes 
it extremely difficult to numerically solve the evolution of the QCD axion from its for- 
mation until the Universe becomes matter-dominated; (ii) in contrast to the ordinary 
curvaton scenario, the total dark matter abundance is fixed by observations. Namely, the 
abundance of the dark matter other than the QCD axion must be adjusted so that the 
total dark matter density satisfies Q c h 2 ~ 0.12. The condition for the adiabatic density 
perturbation, Eq. fl38|) . fixes how to compensate the contribution of the QCD axion. 

Let us recap how the isocurvature density perturbations are generated. The axion ac- 
quires quantum fluctuations of order H in{ /2ii, which results in the fluctuations of the axion 
density. The fluctuations of the e-folding number is through the isocurvature fluctuations 
of the dark matter density. Therefore, in principle, if we know how the axion density 
fluctuates and how the e-folding number fluctuates through the isocurvature fluctuations 
of the dark matter density, we should be able to estimate the isocurvature perturbation 
and its non-Gaussianity. 
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Figure 3: The upper bound on H in{ (solid (green) lines) as a function of the axion decay 
constant f a , for different fractions of the axion to the total dark matter, r = 1,0.1, and 
0.01. We used the constraint from the Planck + WMAP large-scale polarization data. 
The dashed line is for showing the dependence of the constraints on a. 



The homogeneous or spatially averaged dark matter density should satisfy, 



a 



n m + n a ~ o.27, 



(63) 



and it is Q a that acquires spatial fluctuations, while Q m is independent of 5a*. Note 
however that Q m does depend on a*, but its fluctuations do not on 5a*. More precisely 
speaking, the dependence of fl m on a* arises from the requirement that the (spatially 
averaged) total dark matter density reproduce the observed value. However, its fructua- 
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tions is adiabatic and therefore independent of 5a*. Thus, when we expand the e-folding 

number with respect to 5a*, we should use 

dVt c dVt a 
da* da* ' 

keeping in our mind that this is an expansion w.r.t. 5a* and not a*. 

Thus, we may evaluate 5N as follows; 

dN 1 d 2 N 

with 

dN dN dVL a 

da* dQ c da* ' 

d 2 N d^/dnA 2 d^dHla 

~daj ~ afij \da~*~ ) + dQ c ~daJ' [ } 

where we have used Eq. (1M1) . Thus, the power spectrum and the non-Gaussianity are 
given by 

n dN dil a H mi 



" d£l c da* 2tt 

/r> = ±™(™Y + x™y l (!g)(*±Y. ,a 8 ) 



3 dVt 2 \ dVt c J 3 \ dVt c J \ da 2 J \ da* 

The exact solution for the evolution of the Universe which contains radiation and 
non-relativistic matter is parametrized in terms of a parameter £ as 



a, f = 1U +2V1U, (69) 

Hot = ^ c £ 3 + 2 v / ^£ 2 , (70) 

where a s f is the scale factor (it should not be confused with the axion a), t is a cosmological 
time and Q r denotes the density parameter for radiation. In the matter-domination limit, 

3\* „i,„ ,2 



*./(*) - (g) fi "(^)3. (71) 



The uniform-density slicing corresponds to a constant t surface. Thus, 

dN 1 

<Xl c 3f2 c 

d 2 N _ 1 

~d9 2 . ~ ~3fi!' 
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(72) 
(73) 



in the deep matter-dominated era. Finally we obtain useful expressions for Ag and f^ L : 

1 dVt a H mi d\n.Vt a H ini 
\l c da* 2tt ov* 2%j a 

Also) _i I o (® ^ a ' i®® 10 



/NL ^ ^^\dal)\da. 

Thus, given the axion abundance as a function of a*, we can estimate the power spectrum 
A s as well as non-Gaussianity parameter /^ L ' . If Vt a ex a*, these expressions reproduces 
the known results, A s = (2/r)(i7 inf /27ra*), and f£l° ] = -l + l/(2r) (see (J52JI). Note that 
our analytic result for /^ L given in fl6T|) is expressed in terms of the parameters at t — 
tosc, while the above estimate ( 1751) depends on the final axion abundance. The advantage 
of the analytic results in the previous subsection is that one can easily understand the 
behavior of the isocurvature perturbation and its non-Gaussianity in terms of the axion 
dynamics and the potential. On the other hand, the expression ( 175]) is more useful when 
the (approximate) analytic expression for the axion abundance is known. It is also suitable 
for numerical calculations. 

In Fig. HI we compare two expressions for /^ L \ (I6T|) and (175]) . for /„ = 6x 10 9 GeV. 
They agree with each other well for a wide range of the initial condition within about 
10%. We can also see that the expression /^ L = — 1 + l/(2r) is indeed valid only near 
the origin. 

5 Conclusions 

In this paper we have extended the analytic method developed by KKT to a time- 
dependent potential and then applied it to the QCD axion. We have derived the analytical 
expressions for the isocurvature power spectrum (1601) and its non-Gaussianity parameter 
/nl Q6ip . Interestingly, the power spectrum significantly increases in the hilltop limit, 
while the non-Gaussianity increases only mildly. Specifically, /^ L becomes a few tens 
in the hilltop limit when the axion is the dominant component of dark matter. If the 
axion is a subdominant component of dark matter, it increases as l/2r — 1 (see (|62p ). 
The upper bound on the inflation scale H m f becomes extremely tight in the hilltop limit; 
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Figure 4: The solid (red) lines show the semi- analytic result, while the dashed (blue) 
lines represent the numerical results obtained by using ( 1751) . We set f a = 6 x 10 9 GeV. 
Both agree well with each other for a wide range of parameters within about 10%. For 
comparison, the dotted line (yellow) show /^ L = — 1 + l/(2r), valid when the potential 
is approximated with a quadratic potential. 



-f^inf ;$ 1 GeV for f a = 10 10 GeV when the axion is the dominant component of dark 
matter. (See Fig. [3]). Our analytic expressions for the axionic isocurvature perturbations 
will be useful for a probe of the QCD axion dark matter by using the CMB data. 
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A Another derivation of upper bounds on H- m { and 
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NL 



evaluation of k " 



Here let us estimate the isocurvature constraint on H in { by using the analytic expression 
(17411 and a known analytic approximation of the anharmonic effect. To this end, we use 
the axion abundance |19j . 

/ f \ L184 
Q a h 2 ~ 0.195 6?F(6lO( I5i ^) , (76) 



with the anharmonic coefficient |20 



F(9* 



In 



1_£ 

J- _7 



-| 1.184 

(77) 



where we have changed the exponent so as to be consistent with the axion abundance, 
although this change does not affect our results. Substituting these into ( 174|) . we obtain 
the resultant isocurvature perturbation which is constrained by (1561) . We can express the 
constraint as an upper bound on H in { as a function of f a for fixed r (i.e., #* is given as 
a function of f a ). The result is shown in Fig. [5j where the upper bounds are shown for 
r = 1,0.1, and 0.01 from bottom to top. We can see that the constraints are consistent 
with what we obtained by using the semi-analytic formulae (see Fig. [3]). Similarly, we 
can evaluate /^ L by substituting the above expressions into (1751) . See Fig. [61 The 
enhancement toward the lower value of f a is due to the anharmonic effect represented by 

a in am. 
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large-scale polarization data. 
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